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Abstract 



The aim in this note is to define an algorithm to carry out minimal curvature 
spherical harmonics interpolation, which is then used to calculate the Laplacian 
for multi-electrode EEG data analysis. The approach taken is to respect the data. 
That is, we implement a minimal curvature condition for the interpolating sur- 
face subject to the constraints determined from the multi-electrode data. We 
implement this approach using spherical harmonics interpolation. In this elegant 
example we show that minimization requirement and constraints complement each 
other to fix all degrees of freedom automatically, as occurs in gauge theories. That 
is, the constraints are respected, while only the orthogonal subspace minimization 
constraints are enforced. As an example, we discuss the application to interpolate 
control data and calculate the temporal sequence of laplacians from an EEG Mis- 
match Negativity (MMN) experiment (using an implementation of the algorithm 
in IDL). 
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1. Minimization of Curvature with constraints 

This paper focuses on an algorithm for minimal curvature interpolation on the 
surface of a sphere (that is, a 2D problem). The interpolation will then be used for 
a Laplacian calculation. This is straightforward if spherical harmonics are used for 
the interpolation. 

There are many ways to interpolate data on a surface. Here we choose to 
respect the original data points and minimize the curvature of the interpolating 
surface. Note that this may not always be a good idea, as the data may itself be 
"contradictory" due to noise impact. Nonetheless, in the present situation this is 
not an issue. We will be dealing with scalp voltage measurements obtained from 
an even grid of electrodes. No matter how noisy the data is, there always exists an 
interpolation solution for a given set of measurements (e.g., there will not be two 
contradictory measurements at the same electrode site). 

We will thus impose full respect for the original data to fix part of the interpo- 
lation parameter space, and fill the rest using the constraint of minimal curvature. 
The goal is to carry out an interpolation with 

N I 
1=0 m=-l 

minimizing the total curvature: 

N I 

dS |A$| 2 = $([a, m ]) = E E + !) 2 |«U 2 - (1-2) 

1=1 m=—l 

The Laplacian, in spherical harmonics, is given by the simple relation 

N I 

Avty, 9) = -- E E l ( l + !) a i™ Yi, m (d, VO, (1-3) 

r l=0m=-l 

where r is the radius of the spherical head (we used 11.36 cm for an average head). 

Since the array v is real, and since (Yi m )* = Y\ _ m , it follows that (a/ m )* = 
a\- m . This point is important if we want to write the problem using complex 
variables and treating a^ m and (a^ m )* as independent quantities. We have to write 
the problem minimizing something that carries this symmetry forth. Thus, 

X 2 = a) -Ba + {v- Ma) ] ■ X + A f • (v - Ma), (1.4) 

where f denotes hermitean conjugation, is a valid minimizing functional. The fact 
that the interchange of a and a* leaves x invariant means that both will be treated 
equally and extremization with respect to either will yield the same equation, albeit 
in complex conjugate form. 
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The code implementation of this is a bit more complicated for and will forego 
for now the complex approach. In jjj there is an interesting section explaining how 
to translate a complex problem into a real one. 

To translate the problem to the real domain, here we work with the expansion 

N 



1=0 



a\ Y lfi {e,i>)) + a r lm Re[Y ltm (e^)]+ai m Im[Y hm (6,n 



m=l 



(1.5) 



Here v is the array of voltage values at each electrode position. In the present work 
we will have 31 such electrode values. The equation we want to solve (the hard 
constraints) can be written in the form 

v = Ma. (1.6) 

To give a pictorial description of this matrix, M has rows of the form \Yi m (9i,i/Ji)]' 

( Yw(p x ,ik) ii,o(0i,^i)) Re^i^Vi)] ImfOV^Vi)] \ 

Y ,o(e 2 ,ih) Y lt o(02,ih)) Re[Ti,i (03,^2)] M(Fi,i(0 2j ^ 

M 



V ^0,0(631^31) *l,o(031,^3l) Re[Yi,l (031,^31)] Im[(Fi,i(0 3 i^3i)] / 

We will use lower case letters to denote vectors (such as a and A, our unknowns, 
and v, the potential data vector, a 31-vector). Then we have to be extra careful 
with the dimensions of our matrices. 

Let us refer to the original equation M a = v. These are in fact 31 equations 
for an infinite number of unknowns — not a promising perspective. So we change 
strategy! The goal in our game is to minimize the curvature 

C=^a t -Ba, (1.7) 

where B is a symmetric matrix, while respecting the data (the constraints). Note 
that without the constraints in the previous equation, the solution to this mini- 
mization problem is simple: ao is free, the rest must be all zero. We now that 
ao represents the mean voltage, so it is well determined by the constraints. The 
problem we are posing is therefore well defined: there exists a unique solution. 



1.1. The functional and the equations 

Our problem now is to minimize 

x = i a T • Ba + A T • (y - Ma) (1.8) 

with respect to the unknowns a and A. 

These are the players, and we need to present them. First of all, a is the spherical 
harmonics coefficient tensor, the object we are really after. This, in principle, is a 
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oo-array. Then there is v, the aforementioned 31-vector of voltage data values 
measured at a fixed time t. Now, this means M is a 31xoo matrix (the notation is 
" output" x" input") dimensions. Now, let us summarize: 





a 


G 


R°° 




v, A 


G 


R 31 


M 


R°° 




R 31 


M T M 


R°° 




R°° 


MM T 


R 31 




R 31 


B 


R°° 




R°° 



To proceed, let us expand Equation 1.8 



1 i 

X= 2 a 



Ba + A T • - MA) 



-ciiBijCLj + Xi(vi - Mijttj), 



(1.9) 
(1.10) 



and now differentiate to get the stationary point: 

— ^ikOijaj + -aiD ik 



AfAf, 



11.11) 
11.12) 



Differentiating with respect to A yields again Equation |LG|. Our new equations are 
now 



Ba = M T \ 
Ma = v. 



(1.13) 
(1.14) 



Now we have to solve them. The main point is to be careful with the dimensions of 
all the objects — the situation is a bit trickier than usual. 



1.2. Solving the equations 

Let us now add M T times the second equation to the first one (note that, perhaps 
for the first time, we are not adding pears to apples): 

(B + M T M)a = M T (X + v). (1.15) 

This step is important, now we have on the left a matrix of full rank. 

Now, from the second Equation in |1.13| we see (assuming that MM T has full 
rank, as it should and we shall see below) 

A = (MM T )' 1 MBa. (1.16) 
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This we can take and plug into the previous result, 

(B + M T M)a = M T (MM T y 1 MBa + M T v. (1.17) 

Hence 



B + M T M - M T (MM T y 1 MB a = M T b, (1.18) 



rTi 



an important result, and 
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B + M M - M (MM ) MB M 1 b = Q~ xn M b, . (1.19) 
with 

Q nxn =\b + M t M - M T (MM T )' 1 MB] . (1.20) 

Now, can be this be simplified? Conceptually yes. The operator M T (MM T )~ 1 M 
is a projection operator. Let us describe it in more detail because this situation is 
a nice example between constraints and minimization requirements. 

1.3. M and SVD 

To try to see how to simplify, we will use the Singular Value decomposition. Recall 
that a very powerful conceptual and practical tool for studying this problem is 
provided by the Singular Value Decomposition theorem of matrices (SVD)Q, which 
states that given (here m < n, m = 31 and n = oo) any to x n matrix M mxn , there 
exist essentially unique matrices U mxn (m x n), W nxn (n x n), and V nxn (n x n) 
such that 

M mxn = U mxn W nxn V nxn - (1.21) 

These matrices have further properties: W nxn is diagonal, with entries bigger or 
equal to zero, V nxn is orthogonal, 

KxnCn = V T V = l nX n, (1-22) 



and the nonzero columns of U mxn form also an orthogonal matrix(C/ mX n^ r Jxm 
Imxm and 



UnxmUmxn — Pnxn ~ diag(l, 1, .., 1, 0, .., 0) nxn . (1.23) 

There as many nonzero columns as the rank of M as there are nonzero diagonal 
entries in W. 

The power of this decomposition theorem is that it tells us what the kernel and 
range of A are: the kernel of A is spanned by the columns or rows of V which 
correspond to the zero diagonal elements of W, and the range is spanned by the 
columns of U which correspond to the nonzero diagonal elements of W. 
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1.4. MM T 

For instance, consider 

MM T = UWV T VWU T = U 31xoc Wl xoo Ul x31 (1.24) 

Hence, one is tempted to write 

{MM T Y l = U mxn W-? n U? xn , (1.25) 

which is ill defined, since W has nonzero diagonal entries. Let us define, in this 

operation, that 1/0=0. That is, all the diagonal entries of W which are are left 
alone under inversion. Does this work? 

T 



MM T (MM T ) 1 = U mxn W mxn U T nxm U mxn W~$ n U 7 n 



= U mxn Wl xn diag (1, 1, .., 1, 0, .., 0) nxn W n ^ n U, 



nxm 



— U mX n diag(l, 1, .., 1, 0, .., 0) TtX n U nxm 
= ^mxrn (1-26) 

so this is indeed the inverse. 

1.5. M T M 

What about M T Ml This is 

M M = V nxn W nxn U nxm U mxn W nxn V nxn 

= Vn,nW^ xn V n T xn . (1.27) 

This is a large matrix of rank 31. 

1.6. Projection operators: P = M T (MM T ^j 1 M 
Now, 

m t (mm t ) _1 M = V nxn W nxn XJl xm U mxn w~l n ul xm U mxn W nxn V^ xn 

= Vnxn W nxn W nxn Wnxn Vnxn 

= Vnxn diag(l, 1, .., 1, 0, .., 0)nxn V^ xn 

= Pnxn- (1-28) 

This projection operator maps the columns of V corresponding to nonzero entries 
in the diagonal of W to themselves (acts as the identity), and to zero the other 
ones. It projects, therefore, the linear n-space into and m-subspace. 
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1.7. Projection operators: P = 1 - M T (MM T ) 1 M 
Now, recall, 



Qn> 



B + M T M - M 1 (MM 1 )~ L MB 



T\-l 



But 



B - M T (MM T )- L MB = [l nxn - P nxn ]B 



T T\-l 



and 



[1.29) 
[1.30) 

:i.3D 



Q nxn = [M T M + P nxnJ B . 

note how this matrix has full rank, and how the equations have worked out so 
that only the "free" subspace of M T M is affected by the curvature equation. The 
minimizing set of equations are used to complement the constraints to fix a unique, 
well defined solution. 

Let us emphasize that P and P = 1 — P are both projection operators: P 2 = P. 

2. Complementarity: the physical and gauge sub-space 



We can now rewrite equation |1.18| as 

Q-a = b' 

where 



Q = M M + PB, 



(2.1) 
(2.2) 



and 



b' = M T b. (2.3) 
To better interpret these terms, let us work now in the U and V basis. In this basis, 
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(2.4) 



as one would expect. In the second term we schematically show that the projection 
operator "deletes" the impact of B in the physical subspace. What about b'l This 
is similarly given by 

( h 'i \ 



b> 2 

y 



V J 



(2.5) 
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By "physical" subspace we mean the sub- vector space affected by the hard con- 
straints associated, in this example, to the measurements we choose to respect. The 
complement of this sub-space is called, in other contexts, gauge sub-space. It is the 
free-floating part of the solution space, the one that we need to somehow fix. 

Thus, we see that the method of minimization subject to constraints leads to 
a nice interpretation in which the physical degrees of freedom (associated to the 
constraints) are not affected by the minimization extra requirement, while the gauge 
degrees of freedom are fixed by the curvature minimization requirement. In the 
language of fiber bundles, which is appropriate for handling gauge theories and 
which we could have used in the above discussion, by the choice of B we are choosing 
a particular (gauge) fiber bundle section. 

This method implies that we trust fully the constraints, that is, we trust fully 
the measurements, as if there was no noisy. This is, in general, not a good strategy, 
as the data may be very noisy and inconsistent. For this purpose, following our 
intuition as based on the previous discussion, we can generalize a bit Q to write 

Q = M T M + PB + vPC, (2.6) 

where v is a "tuning" parameter and C is a new (normalized) constraint affecting 
now only the physical subspace (note that tuning the middle term has zero impact) . 
An immediate choice is C = B, of course. 

The first term in Q is acting only in P-space. This encodes the limited in- 
formation available from the data. The second one encodes the information from 
curvature minimization requirement as projected to P-space. The third, new, term, 
encodes any additional information we may want to add to "smooth" or regularize 
the information in the first term, which may be noise or otherwise unreliable. 

How much should be added to the physical part? As much as needed to fix the 
solution, but no more. That is, if we have some a priori knowledge on the noise 
characteristics associated to the P-space constraints, we can evaluate the entropy 
associated to this condition. We then need to add enough information to bring the 
entropy down to the desired (or needed) value. 



3. The algorithm 

The goal is to implement code to solve 



a = 



with 



B + M T M - M t (MM t )- 1 Mb] M T b = Q~l n M T b, (3.1) 



Q nxn =\b + M t M - M t {MM t )- 1 Mb \ . (3.2) 



Once the coefficients are available, a new M is constructed from a interpolated set 
of electrode positions, iM. This matrix is specified here by a set of 40x20 positions 
(the larger number for longitude positions). 
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Using the interpolating matrix, it is then simple to compute the interpolated 
potential field and its Laplacian. 



4. Analysis of sample experimental data 

Here we give a simple example using an implementation of this algorithm devel- 
oped at Starlab. A MMN sample data set from an Experiment carried out at the 
University of Barcelona Neurodynamics Laboratory is analyzed for illustration. For 
the purposes of the present discussion it suffices to mention that the 32 electrode 
(plus a reference) data set corresponds to an average of MMN EEG from 17 control 
group subjects, and that it corresponds to 100 ms prior to an auditive stimulus up 
to 500 ms later. For more information see and forthcoming publications. Here 
we show results for lmax=20. 

4.1. Control group 

This is the output for the data interpolation (32 electrodes, time in ms): 




time= 100 time= 110 timc= 120 timc= 130 



^ I ex 20 1 I n 20 j ! cn 20 j 

^^^^T - • r; ■ Wgf - o - -io ^M^^r 

^^^^^ u _9n ^^^^^ u _9n ^^^^^ u _9n ^^^^^ 



-20-10 10 20 
time- 1-10 




-20-10 10 20 
time= 180 



20 r 

10 



-10 

-20 L 



-20-10 10 20 
time= 150 



20 f 
10 


-10 
-20 L 



-10 10 20 
time= 190 



-20-10 10 20 
time= 160 



-20-10 10 20 
time= 200 



-20-10 10 20 
time= 170 



-20-10 10 20 
time= 210 



-20-10 10 20 
;irre = 220 



-10 10 20 
time = 230 



1-10 10 20 
time= 240 



-20-10 10 20 
time= 250 



And this is the resulting Laplacian: 
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-20-10 10 20 
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-20-10 10 20 
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-20-10 10 20 
time= 1 10 




-20-10 10 20 
time= 150 




-20-10 10 20 
time= 120 




-20-10 10 20 
time= 160 




-20-10 10 20 
time= 170 



20 
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-10 
-20 




-20-10 10 20 
time = 180 



o 20 
- 10 



, ; ~ ^ 20 _ , _ Ln 20 

^fc^MF ■ 3-10 ■ 5-io N^^Br 

— I " -201 w ! J -201 ^^^^ 



-20-10 10 20 
time= 190 



-20-10 10 20 
time= 200 



-20-10 10 20 
time= 210 



20 ' 
1 ■ 



-10 ' 

-20 . 



-20-10 10 20 
:in^e= 220 



-20-10 10 20 
time= 230 



-20-10 10 20 
time= 240 



ID 10 



-20-10 10 20 
time= 250 



4.2. Curvature index 

This a graph of the curvature of the fitted minimal surface, an interesting measure 
(the blue curve, the other curve refers to another experimental group). 
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5. Analogy to Constrained Mechanics and Gauge Theories 

In this set of notes I intend to introduce very quickly to the Lagrangian for- 
malism in physics, the impact of symmetries of the action to the solution space, 
their relation to gauge theory, and how all these strange things happen to relate 
to the interpolation problem (for more on this and further references, see xxx-gr- 
qc9806058). What all these have in common is their origin in a minimization (or 
extremization, to be precise) problem. 

5.1. The Lagrangian in physics and x 2 

It is interesting to form an analogy to the theory constrained theories in classical 
mechanics. 

Given a system with ndegrees of freedom, each specified by a coordinate <?;(£), 
we can obtain its dynamics by minimizing the functional (called the action) 

S[q t (t)} = f 2 L( qi (t), qi (t))dt, (5.1) 

that is, one looks for the set of qi(t)such that given that at t± and i2the values of the 
coordinates are to be held fixed, one varies the functions so as to find the minimum 
of S. In general, and although this may seem less familiar to you, one can rewrite 
this as 

SMt)] = f 2 /" 2 q t (t)M^(t,t%(t') dtdt' - f 2 V( qi (t)) dt (5.2) 

If the last term is a quadratic (as it happens in the case of doable physics, e.g., the 
harmonic oscillator) , we find that we are trying to minimize a quadratic functional, 
just as was the case in the previous sections, where the goal was to minimize the 
curvature. There is no single surface of minimal curvature, as a constant term can 
always be added, for example. 

Allowing for linear terms, then, we are trying to minimize something that looks 

like 

S = q T Mq + q-j + c (5.3) 

where q here stands for a long vector in which there is an entry for each time 
parameter (think of time as a discrete index if this helps)and for each possible value 
of the index i. By a symmetry we mean a transformation of the quantities q — > q'so 
that S — ► S: i.e., so that the action is invariant. For instance, and to simplify 
things, imagine that the action we want to minimize is 

X 2 {x) = {y-Axf (5.4) 

where x, which we will call here the state vector or state, for short, is an n- vector, 
and y is an m- vector (this equation may look familiar to you!). Let us now use the 
Singular Value Decomposition theorem to rewrite 

A = UWV T (5.5) 
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Recall that W is a n x n diagonal matrix with entries bigger or equal to zero, V is 
an orthogonal n x n-matrix (VV = T VV T = I n xn), and that U is and m x n-matrix 
with orthogonal columns (U T U = I n xn)- 

Imagine now that the last m diagonal entries of W are zero. What does this 
mean? It means that the two states x and x + ajV^are equivalent as far as \ 2 ls 
concerned, where by Vj we mean the j-th column of V: 

X 2 {x) = X 2 {x + aiV l ), i = l,...,m (5.6) 

The fact that the m coefficients aj are all independent of each other is what makes 
this transformation symmetry a local symmetry in the physics lingo. It tells us that 
of the original n degrees of freedom, only n — m really matter. These are called 
the physical degrees of freedom. The rest are called gauge or unphysical, degrees of 
freedom. These are the source of our inversion problems: the minimization problem 
does not have a unique solution. In physics one fixes this in some way or another, 
usually by adding a so-called gauge-fixing conditions. This means to pick a set 
of coefficients ctj, and this is usually done by adding a set of m equations to the 
problem that fix these. Let B be a rank-m n x n-matrix. Then we'd like to add 
such a condition to the problem. Notice the condition Bx = is to fix the gauge 
degrees of freedom then B has to have maximal rank in the gauge subspace (i.e., it 
must completely span the gauge subspace) . 

If we believe in the data, one must avoid increasing x 2 while fixing the gauge 
degrees of freedom. It is possible, as we now discuss, to fix these in a reasonable 
manner, without getting away from the data. We have in fact discovered how to do 
that in the previous discussion, but it may be useful to repeat the exercise in this 
context. 

Suppose then that for some physical reason we want to ask that Cx = in 
order to fix the gauge degrees of freedom, but that we have been a little naive and 
have not worried about the rank of C. Proceeding as is usual, then, we choose to 
minimize the functional 

X 2 (x) = (y-Ax) 2 + X(Cx) 2 . (5.7) 

Upon minimization we obtain the equation 

(A T A + \C T C)x = A T y. (5.8) 

Let x = Vx' . This is the state in "SVD coordinates". The first n — m entries in 
the vector x' are physical coordinates, the rest are gauge. Then, using the fact that 
A T A = VW 2 V T , we can rewrite this equation as 

V T (A T A + \C T C)VV T x = V T A T y, (5.9) 

or 

(W 2 + \V T C T CV)x' = WU T y. (5.10) 
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This equation reflects the fact that without the additional constraint (i.e., set A = 0), 
the physical degrees of freedom are well-fixed, while the gauge ones are not. If we 
want to avoid disturbing the physical part of the state, we now project the constraint 
into the gauge subspace: V C T CV — > P gauge V C T CV, where Pgauge is a diagonal 
matrix with zeros everywhere except the last m last diagonal elements. Now the 
solution to 

(W 2 + \P gauge V T C T CV)x' = WU T y (5.11) 

is as before as far as the physical degrees of freedom are concerned. If the rank of 
PgaugeV C CV is m we will have fixed a complete solution. To see more clearly 
what has happened, let us factorize x' = x p h y + x gauge , where x gauge ~ oiiV % is the 
projection of x into the gauge subspace , we can rewrite the above equations as 

W 2 x phys = WU T y, (5.12) 

which fixes the physical degrees of freedom , as it should, and 

PgaugeV C CV {x p hy + X gauge ) =0. (5.13) 

In this way we ensure that we are not imposing any constraint equations beyond 
those that we are really allowed to. The matrix P gauge V T C T CV is to have rank m. 
If C is fully ranked, then it will. If not, we may still be ok. We can now go back to 
the original coordinates and write (notice that A is irrelevant now) 

(A T A + VP gauge V T C T C)x = A T y (5.14) 

Another way to reason this result is the following. As we mentioned, the physical 
degrees of freedom, P p h ys x, are fixed. We want to really find now the set of gauge 
degrees of freedom such that (C = V T CV is the expression of the constraint in the 
SVD coordinates) 

(C {PgaugeX + Pphys x )) (5.15) 

is minimum. If we vary the gauge degrees of freedom in this functional and set it 
equal to zero we find the equation 

(^PgaugeC C Pgauge PgaugeC C Pphysj X = 0, (5.16) 

or PgaugeC 17 ' C x' = 0, since Pgauge + Pphys = 4xw In the original coordinates, this 
reads V P gaU geV T C T C x = 0, as we wrote before. 

By construction, then, there exists now a unique solution to the problem 

(A T A + VP g augeV T C T C) x = A T y. (5.17) 
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